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ABSTRACT 

This paper presents a formulation for identification of 
linear multivariable systems from a single set of input- 
output data. The identification method is formulated with 
the mathematical framework of learning identification, 
by extension of the repetition domain concept to include 
shifting time intervals. This approach contrasts with 
existing learning approaches that require data from 
multiple experiments. In this method, the system input- 
output relationship is expressed in terms of an observer, 
which is made asymptotically stable by an embedded 
real eigenvalue assignment procedure. Through this 
relationship, the Markov parameters of the observer are 
identified. The Markov parameters of the actual system 
are recovered from those of the observer, and then used 
to obtain a state space model of the system by standard 
realization techniques. The basic mathematical 
formulation is derived, and numerical examples 
presented to illustrate the proposed method. 


INTRODUCTION 

The aim of learning system identification is to 
provide methods to improve identification of system 
parameters as new data in the form of inpul-output 
measurements are available. New information regarding 
the system characteristics may come from multiple 
experiments. For system identification of flexible 
structures, multiple experiments arc usually performed to 
develop or improve a mathematical representation. Due 
to structural complexity and data irregularities such as 
slight non-linearities, instrumentation errors, background 
noises, and repetitive disturbances, multiple tests are used 
to reduce the irregularity effects on the identified model 
parameters. The conventional approach is to average data 
sets from multiple experiments with the hope that the 
averaged data will reduce the irregularity effects. In 
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learning identification, new information from successive 
experiments arc used to effectively improve current 
identification result. In fact, current work in learning 
identification falls within this conceptual framework. 
They require the availability of different input-output sets 
of data from different experiments of generally short 
duration. However, new information about the system 
need not come from new experiments, but rather it may 
be derived from a single experiment of extended duration. 
This motivates the development of an identification 
algorithm using a single set of input-output data. 
Originally motivated by the learning problem, the 
algorithm is derived using the mathematical framework 
and techniques of learning control and system 
identification. Refs. [1-9], and related in general concepts 
to identification methods proposed in Refs. [10,11]. In 
particular, this requires the extension of the concept of the 
repetition domain to shifting time intervals. 

In usual identification techniques, the system time 
domain parameters are determined from input-output data. 
In learning system identification, however, the parameters 
of interest to be identified are the Markov parameters. 
Once this step is completed, standard realization 
procedures can then be used to realize the system time 
domain parameters. This shifting emphasis on 
identification of the Markov parameters offers some 
advantages. First, there is no ambiguity in the dimensions 
the Markov parameters. Second, in the learning 
formulation, the Markov parameters arc related to input- 
output data by a simple linear relation, hence many 
existing techniques can be applied. Third, for a given 
linear system, the Markov parameters arc unique and 
invariant with respect to any coordinate transformation of 
the state vectors. 

In this paper, a treatment of this problem is presented 
to identify a linear multivariable system in state space 
format by first identifying its Markov parameters. From a 
single set of input-output data, direct solution of the 
system stale space matrices is non-trivial for a general 
system of completely unknown characters. However, 
when additional information about the system is imposed, 
the mathematical problem becomes simplified. For clarity 
of exposition, identification procedures arc presented for 
the following three cases of increasing complexities. 
First, the order of the system is known, and is equal to 
the number of outputs. Second, the order of the system is 


1 


not known, but the system is known to be asymptotically 
stable. Third, the order is known, but no assumption on 
the stability of the system is made. The resultant scheme 
is recursive coupled with an eigenvalue assignment 
procedure, and is based on techniques developed in the first 
two cases. The algorithm has an embedded observer 
structure with pole placement. It is emphasized here that 
the role of the observer structure is not to provide 
estimates of the system states for identification, but rather 
to provide by design a set of asymptotically stable auto- 
regressive moving average equations whose parameters 
can be identified. These parameters contain in them the 
desired information about the actual system. The initial 
assumption regarding the system order can be later 
removed by an iterating process, or by knowledge of an 
upper bound on the effective order of the system. 

The Markov parameters of a linear system in state 
space format arc related to the system impulse response 
functions which can be used in modal identification. From 
the identified Maikov parameters, the modal parameters of 
the system such as natural frequencies, modal damping, 
and mode shapes can be deduced by standard realization 
procedures, e.g., the Eigcnsystem Realization Algorithm 
(ERA), Refs. [12,13], 

STATEMENT OF THE PROBLEM 

Consider a general discrete multivariable linear 
system expressed in state space format as 


i-0+ 1) 

x(t) = A ip x(p) + X A iz ^ X) Bu(x+p) (3) 

t*0 

Introduce a repetition variable j , and a new time step 
variable k , so that the general solution to Eq. (1) for the 
first, second, and all subsequent intervals can be expressed 
as 

<*-(/>♦ i) 

x(i) = A iip x(jp) + X A‘ T<Jptl) Bu(x+jp) (4) 

T*0 

for / = jp + k, j = 0, 1 , 2, ..., k = 1,2, ..., p. At any time 
step i = jp + k , the state vector x (i) is written as 

x(jp+k) = Xj{k) and xfjp) = x/0) (5) 

and similarly for y(/p+k). Associated with each state 
vector x(0 is an input vector u(/-l) which is also rewritten 
in terms of j and k as 

«(M) = «/k- 1) (6) 

Applying the definitions in (5) and (6) to (4), one obtains 
the following description of the system in the repetition 
domain 

tbi =Axj(fl) + / > k,+i ( 7 ) 


x(k+\) = Ax(k) + Bu(k) 
y(k) = Cx(k) 

where x e R n t y € R q , ue R m . The number of inputs 
m, and of outputs q are known, the order of the system n 
is in general assumed not known, and neither are the 
system matrices A t B t C . Starting from some arbitrary, 
possibly unknown initial state x(0), the system (1) yields 
a sequence of outputs y(i) t when driven by some known 
sequence of inputs #(M), i = 1, 2, 3, ... for an extended 
amount of time. The main objective of the problem is to 
recover the set of Markov parameters of the system 

CB y CAB , CA P l B for a given value of p. For 
purpose of identification, knowledge of a sufficient 
number of the system Markov parameters is adequate to 
deduce a state variable description of the system, hence 
completely characterizes the system of interest 

MATHEMATICAL FORMULATION 


where 

0 ( 1 ) x t (2) ... x T (p)f 

« = O(0) m t (1) ... u T (p- l)f 



P = 


B 

AB 


O 


B 


A pX B A P ' 2 B 


B 


using y/k) = Cxj{k) y the equivalent output description to 
(7) is 


First, the set of known input-output data is divided 
into intervals of p time steps each. For j = 1,2, ..., p , the 
solution to (1) is 

x([) = A i x(Q) + Y J A ixX Bu(x) (2) 

f-0 


yj+ 1 = GAxj+ i (0) + Po uj+ 1 (8) 

where 

X = [y T (\)y T (2)...y T (p)] T 

£ = diag[C C ... C] Po = £P 


For the next interval, / = p+ 1, p+ 2, .... 2 p, the solution 
is written as 
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M*)=|W -l)«./,l P(k - ••• P(*-l)n+/^i*] r 


B o = 


CB 

CAB 


o 


CB 


L CA' 'B CA p ‘ 2 B 


CB 


Note that in the repetition domain, the unknown sequence 
of Markov parameters to be determined, CB, CAB, CA 2 B, 
... appears naturally. One can make the following general 
observations before proceeding to the respective 
algorithms. 

Equation (8), in a succinct form, relates the known 
input-output data in terms of the system Markov 
parameters, and necessarily the quantities CAjc /+ i( 0). If the 

states x>+i(0) are known, then given a sufficiently long 
sequence of input-output data, subject to some appropriate 
conditions on the input sequence, the system Markov 
parameters can be uniquely recovered. In general, however, 
the states xy + i(0) arc not known. Equation (8) when 
written for all available input-output data, represents a set 
of under-determined equations, with the quantities Xj* i(0) 
as additional unknowns. Therefore, the system Markov 
parameters cannot be uniquely determined without 
imposing additional constraints to the set of equations, 
i.c., without assuming additional knowledge about the 
system to be identified. In particular, we consider the 
following cases. 

If the order of the system is equal to the number of 
outputs, and the system is observable, i.c., if C is square 

and full rank, then with Ao=CAC’ 1 , the output 
description (8) can be rewritten as 


where the subscripts on a and p denote the position of the 
elements in Aoand Bo respectively, l = 1,2, .... q,k= 1, 
2,..., p. The /nk-dimcnsional input vector is defined as 

a«,/+t(*-l) = [«i.;>i(0) ... u m ,j+ 1(0) ... u m ,j+\(k-l)] T 

Making use of the above definitions, Eq. (9) can be 
expressed as 


y, .,♦,(*) = aJ(k)yM 0) + p!(k)Wn,Mk-l) (10) 

for all (l,k) pairs. Equation (10) represents a set of single- 
output, multiple-input models in the repetition domain. 
Grouping the unknown parameters together, Eq. (10) can 
be rewritten as 

y/.M(*)= P/ + (*) r j4/*i(*-l) OD 

where 

p!(k) = [a[(k) p?(k)J 

iC^i(*-l>=[^i(0) «£.;♦ i(*-l)] 

The parameter estimates of Pi(k) at repetition j, denoted 

A + 

by Pijih) can be updated recursively by various methods, 
for example, 

i) The Projection Algorithm: 


yj* i = Aoy>+i (0) + Pom+\ (9) 

From Eq. (9) it is clear that subject to some usual 
condition on the richness of input sequence, from a given 
set of input-output data of sufficiently long duration, the 
Markov parameters can be uniquely determined. 

There arc various methods that can be employed to 
solve for the Markov parameters from (9). Here, we arc 
particularly concerned with recursive algorithms. 
Recursive algorithms offer some fundamental advantages 
for this particular problem. Namely, they arc efficient in 
processing a large amount of data on-line. Furthermore, 
recursive methods arc essentially approximation methods 
to solving the problem. As such they can be tailored to 
solve for certain aspect of the problem without solving 
the full problem as would be required by any exact 
method, which is some cases, may require exact solution 
to a difficult non-linear problem. To solve (9) recursively, 
first let the rows of Ao, Bo be defined as 

<z/(k)=[a<*-iV, + ;.i a(t-i)»+j,2 ••• a(* i )<.+/,» J 7 


pt/k) = ph-\(k) + ajtC,,j{k- 1) 


yjJthjCfa QPO-i(*) 

. 1 + tU T j(k-\)tUj{k-\) 


for all (/, k) pairs, 0 < a j < 2. 
ii) The Least Squares Algorithm: 


pij(k) = ptj-\(k) + Rj 2{k)iC n , j{k - 1 )Apij(k) 


Apijk) = 


yi.j(k) - lU T i<k-\)pijA{k) 


Ll + i£ T j<k-\)Rj-2(k)Cj(k-l)\ 


Rj ■(*) = 


RMk) - 


R^Jc^Jk- l)i£, r /M)/?;.2(fc) 
. 1 + lU T j{k-\)Rj-2(Jc)C,j{k-\) 


R.\(k) = almkx mit, for some large a>0. In actual 
implementation of this algorithm, covariance resetting 
may be employed to speed up the convergence of the 
parameter estimates. 
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Other recursive algorithms can be easily incorporated 
within this framework, and arc not considered here for the 
moment. One can analyze the above algorithms to 
establish the necessary condition on the input sequence to 
guarantee that the estimated parameters do indeed converge 
to the true values. For the present development, another 
class of problems where recursive approximation methods 
can be used to recover the system Markov parameters is 
considered. 


0 CA^B ... CA 2 B CAB " 

0 0 CA^'B ... CA 2 B 

_0 0 ... 0 0 

after neglecting all terms involving p or higher powers of 
A . Then 


Consider a special class of problems, where the 
system matrix A is known to be asymptotically stable. In 
particular, p can be chosen to be sufficiently large such 

that A k can be approximated to be zero for k > p. Then 
the sequence of Markov parameters CB f CAB , CA P B , ..., 
CA p] B can be recovered by the following method. First, 
Eq. (8) is rewritten here for clarity and convenience 

Xf+t = CAxj + 1 (0) + Pouj+\ (12) 

We seek to write an expression for*y +1 (0) from (7). Note 
that (7) can be rewritten as 

% = &Xj.x(Q) + Puj ( 13 ) 

which yields 

Xjip) = A p Xj.\ (0) + \^ p X B ... AB b}% (**) 
Defining 

P(p) = [a p 'b A p ' 2 B ... AB fl] 

and noting that xj(p) = ^ +1 (0), Eq. (14) becomes 

x>+,(0) = A p xj (0) + P(p)Uj = P(p)m (1 5) 

after making the approximation A p - 0. Substituting (15) 
into (12) yields 


[CAP(p) />„] = 


0 CA pl 

B - CA 2 B CAB 

CB 

0 

0 

CA pl B ... CA 2 B 

CAB 

0 


: 

: 

: CB 0 

0 

0 CA pA B 

CA p ' 2 B 

... CB 0 

- 

0 

CA pl B 

CB 

Defining 


Po(p) = [cA p] B ... 

ca 2 b 

Cfi] < 18 > 


to be the matrix of Markov parameters to be identified, 
Eq. ( 1 6) thus becomes 

yj+ i(*) = Po(p)uJ+\(k-\) > *=1.2 p (19) 

where i(*-l) is a mp-dimcnsional input vector defined 
as 


[uf(k) «;,(0) «;,(d ... «j,(*.i)f 

To estimate P 0 (p) recursively, the rows of P Q (p) can be 
recursively updated in parallel. Let y£k) denote the /- th 
output at lime step k , and pi denote the column vector 
formed by the /-th row of P Q (p ), Eq. (19) becomes 


*+t = £AP(p)iij + PoUj+x 


-[CAP ip) 



(16) 


Consider ihc product CAP(p) 


CAPip) = 


CA 




CA 2 
CA I'” 1 

[a^B A p ' 2 B 

... AB 

B] 07) 

CA P 




CA P B 

ca^'b - 

CA 2 B 

cab' 

CA^B CA"B - 

CA'B 

ca 2 b 

CA^ 'B CA ^ ' 2 B ... 

CA^'B 

CA p B, 


yi.i+m = pJ^{k-\) ( 21 ) 


for / = 1,2, q. Each /-th row of P 0 (p) is now updated 
from input-output data associated with the /-th output at 
lime step k t k = 1, 2, of repetition y+ 1 , y+2, ... . 
The resultant identified parameters are smoothed both in 
time and in repetition. Hence, the identification scheme 
makes use of all available input-output data to arrive at 
one single set of parameter estimates. Once the 
identification is put into the form as in Eq. (21), a 
recursive solution is simple, for example, the recursive 
least squares estimate is 


P*j(k) = p*j-i(k) + Rj-2(k)uf (k-\) Ap*j(k) 


Aplik) = 


yUk)-u; T {k-\)pU(k) 
l + U* T (kA)Rj. 2 (k)^(k-\) 
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*/-!(*) = 


*j-2(k) - 


i +^ T (k-i)Rj.2(k)u;(k-\) 


xtf+l) = Ax(i) + Bv(i) 
y(i)=Cx([) 


Note that in the above algorithm, the identified system is 
required to be asymptotically stable. In the following 
development, this restriction will be removed. 

Consider the system in Eq. (1). It has an observer of 
the form 

*(/+!) - Mi) + *11(0 - Af [y(0 - m 
= (A+MC)x(i) + Bu(i) - My{[) 

Defining the state estimation error x{i) = xfj) - x(i ), the 
equation that governs x(i) is 

*(/+!) = (A + AfC)I(i) (23) 

If A/ is chosen such that A +Af C is asymptotically 
stable, then the estimated state will converge to the true 
state as i tends to infinity. If the system (1) is observable, 
then Af may be chosen so as to place the eigenvalues of 
A +MC in any desired (symmetric) configuration. The 
above is a Luenberger observer, a well-known result. 
However, in our problem, the system matrices A, B, C, 
and hence Af are not known. Multiplying C to both sides 
of (22) yields 

y(i+ 1) = CAx{i) + CBu(i) - CMy(i) (24) 
where for simplicity, 

A = A+AfC 


which when compared to (1), A f B , v(i) play the same 
roles as A, B, u(i) t respectively. 

It is noted here that since Eq. (26) is derived from 
(22), hence it can be interpreted in terms of an observer 
equation. But in fact it is an exact relation which always 
holds true for any matrix Af . Equation (26) can be derived 
by simply adding and subtracting the product My(i) to the 
right hand side of Eq. (1). Recall that if the system (1) is 
observable, then there exists at least a matrix Af such that 
the eigenvalues of A can be assigned in any arbitrary 
symmetric configuration. Such an Af is not necessarily 
unique, but this poses no restriction to the present 
problem. In the following development, the additional 
freedom introduced by the matrix Af will now be used to 
derive a recursive algorithm to identify the Markov 
parameters of the system in (26), i.e., 

CB t CAB , ..., CA^ **, and at the same lime place the 
eigenvalues of A in desired asymptotically stable 
locations. Note, however, that the identified parameters 
are the observer system Markov parameters, but they will 
be used to recover the desired actual system Markov 
parameters. This will be done in later sections. 

Applying the technique developed before for 
asymptotically stable system_to the system in (26), 
assuming for the moment that A is asymptotically stable, 
yields the corresponding version of (19) 

Wj^i(*+1). £ = 1,2, ..., p (27) 


If A is asymptotically stable, then for large i, x(i) tends to 
x(i), hence y(i) tends to y(i), then Eq. (24) becomes 


y(/+l) = CMO + CBu(i) - CMy(t) 
= CAx{i) + C*v(0 


(25) 


where 


b = [b -a/] 


and 


v(/) = 


u(i) 

b(0i 


Poip) and arc defined as 


B °(P) = [cA P ' i B ■■■ CAB 
and 

00 

l£i (*-0 = 

(29) 

[v/(A) ...v,Vl), «ai(0) v;,(d ••• 

^,(M)f 

where 



If A is asymptotically stable, then by the method 
developed above for asymptotically stable system, the 
Markov parameters of the observer equation, 

CB y CAB , ..., CA pA B can be identified. This can be 
easily seen by noting that from (22) and (25) for 
asymptotically stable A and large i, the system in (22) 
becomes 


v /(*> =[«/(*) y/w] 

The eigenvalue assignment procedure can be derived, first 
by noting that for desired (real) eigenvalues of A, we have 
for some T 

A = T a AT (31) 


5 


where 


A, 


A = 


A 2 


0 


0 

A, 


(32) 


Let A..A 2 A n denote the desired eigenvalues of 4 to 

be selected a priori. Then the product CA* B becomes 

cT*fl = crVra p 3) 

= cVfl* 


If the elements of C* and 6" arc written explicitly as 




c n c i2 

— • 
a 

1 




c* = 

• * 

C21 ^22 

•• cl, 






* • 




bn 

hi 2 •• 

b\m 


• •«*, 

fl* = 

b 21 

bn 

.. j,. 

“^21 -^22 *• 



1 

r. • 

b\t ... 

bnm 

2 •• 



then in general, the /-th row of the matrix product in (33), 
denoted Pi , is equal to 


p\ k) = f X ^ c i‘bn . X ^ - X ^ c ‘‘b 


•*-v 

stm * 


*=I /-I i-1 

- X C « m <* * - X *** C l‘ m ‘2 . • • • . - X ] 

1=1 .= 1 <=1 


Let the /-th row of /’„(/>) of equation (27) be denoted as 

pi , then the /-th output at time step k of repetition j+ 1 , 
denoted by y (</+1 (*), can be expressed as 


yi.iMk) = pl y.^i(k-]) (34) 


where 




1) -(/>- 2) 
Pi 


-d) 


P/ 


n (0) l 

P/ J 


Note that each pj k) is a row vector, each with m+<? sums. 
Each sum has n terms. Together each p$ k) has n(m+q) 

terms of the products c*fe/*and c*m* k . First, write for m 
inputs and q outputs 


«/(*) = [«!./<*) «2./*) ••• u m ,j(k )] 
and 

?/(*)= [yi./*) n.){k) ••• y,./*)] 

Note that 


hi./*) x = 
1-1 


[cnhii Cnbi\ - ■ cJ» H \\ 


Af ' 


Hi./*) 


A? 


-l 


With the following simplifying definitions of the 
unknown parameters 

atr = [ C|ji»lr C/jhJr 
fib ~ [cjifflj, C/l/Wj, ... Cijrims 
r= 1 , 2 ,..., m ; s= 1,2,..., 9 
and a vector of the assigned eigenvalues 

a;-' ... w'J (36) 

Equation (34) can then can be expressed as 


y/.y+i(*) = 


«n 

I A B ~'"VX«> + I 


L T =* 

T-0 


p-i 

*-l 

+ cr£ 

X * (i 

"'V/ T)+ X 


T * * 

f-0 


M 


+ «lm 

X * (t 

t+pl) « m .XT) + x i M W.w 


T — * 

T-0 


p-1 

*1 1 

- fin 

X a^Vxd+ X 


r- k 

j 


p - 1 


- fil 

X^‘ 


f — k 

T- 0 


iv« 

1 

- fii ; 

X A (it+p - , W/T)+ X 

T«* T-0 
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or 


Rj-i(k) = RMk) - 


Rj-mW-Vifik-Wrfk) 

. 1 + rl(k-\)Rj. 2 (k)m-i) . 


m 

x A (i -"' V.XT)+ x * <tT V»(*) 

f-i 

T- k T-0 

1 

/>'! *1 

■ I A 

X A (t -^ ,) y«.XT)+ X A (i Tl k>.,(T) 

{■1 

T— k T-0 


This equation can be simplified further by making the 
following definitions 

0r./ + l(*-l) = X A^' V/T) + X ^Vi*l(T) 
T« * T - 0 

V^.y + iOt-l> = X A ( * r+#,1) yx.XT) + X 

T»* T-0 

Then it becomes 

yij*i(k) = ah <tn,j*\(k-\) + • ■ • + aj m <t> m .j*\(k- 1) 

- Pi l ■ ••• ■ Pi<i Yi.j+i(k- 1) 


The above algorithm in (39) identifies the observer 

parameters yi which consists of a n , aj m , Pn< 

Pi 2 Pi r Each a, u a l2 , Pn . Pa . in turn contains 

certain products of the elements of C* and B’, e.g., 

c,*/>u, c,*/n fi These products together with the 

assigned eigenvalues Aj, A z A„ can be used to 

construct the /-th row of the observer Markov parameters 

CA k B. Thus any row of the observer Markov parameter 
for any given k can be computed. Hence a complete set of 

— — P"1 — 

the parameters CB , CAB , ..., CA B can be identified. 
Furthermore, note that Eq. (38) is true for all k = 0, 1, 2, 
.... p ; and j = 1, 2, 3, ... hence the identified y h which is 
then used to recover the /-th row of CA k B is smoothed 
over all k and j. Note that Eq. (27) is simply an ARMA 
model of the system, and the observer Markov parameters 
are precisely the matrix coefficients of this ARMA model. 
Through the embedded eigenvalue assignment procedure, 
the ARMA model is made asymptotically stable by 
design. It is the stabilizing property of this model that 
allows a straight forward application of existing parameter 
estimation techniques to identify its matrix coefficients. 
What remains to be done is to recover the actual Markov 
parameters of the open loop system from these identified 
matrix coefficients. 


i(*-l) 


•[ 


ah 




0«.y+i(A-l) 
VUj+\(k- 1) 


Lv'v.Mtf-Dj 


In the following a procedure to recover the actual 
system Markov parameters is presented for the general 
case of mujtiple-input, multiple-output systems. First, 
recall that A = A +MC, B = [B, -M ], and the Markov 
parameters of the observer system are now known from 
the identification algorithm developed above. For ease of 
presentation, the following definitions arc made 


By making obvious definitions for y J, and 7}+i(A-l), the 
above equation becomes simply 

yi.Mk) = yj />i(*-D (38) 


for / = 1, 2, ..., q. The above set of equations is in linear 
form, with unknown (time-invariant) observer parameter 

vector yi and known "input" vector 7>+i(Jk-l), therefore yi 
can be solved for in one step, or recursively. Any 
appropriate method to solve linear equations can be used. 
For example, the recursive least squares solution to (38) 
can be written down immediately for the /-th row of 

C*A k B *, which corresponds to the /-th output, as 


yi.jik) = ytj- \(k) + Rj- 2 (k)rj(k-\)Anj(k) 


Ayt,j{k) = 


y,.j(k) - rf(k- ljg^iW 

1 + rl(k-\)Rj-2(k)rj(k-\)_ 


(39) 


Y 0 = CB = [ cb> .cm] 

"Iff’, yH 

(40) 

Fi B ] 

and similarly for Y 2 , T 3 , ... From the first equation in 
(40), the first Markov parameter of the actual system, 
To = CB t is simply 


r 9 -n" 


(41) 


Next, consider the product C A B 

CAB = C(A + MQ[B , -M ] 

= [CAB + CMCB , -(CAM +CMCM )] 

-to 1 '. ?H 
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Hence, ihe second system Markov parameter Y\ = CAB 
can be computed as 


Y x = CAB 


= rt i ) . n 2) Y o 


(42) 


2 

To obtain an expression for Y 2 = CA B, consider the 
product Y 2 = CA^B 

CA 2 B = C(A+MC) 2 [B ,-m] 

-In" , n n ] 

Thus 

Pz 0 = C(A + MC) 2 B 

= C A 2 B + (CM)CAB + ( CAM + CMCM)CB 

Hence 


Y 2 = CA ‘ 


= y 2 1) + Yo 2) Yi + yI 2) y 0 


Similarly , an expression for Y 3 can be obtained from 

ac 


as 


Y 3 = CA’B 

= y] 1) + y1 2) y 2 + Ff Vi + y1 2) Yo 


( 44 ) 


By induction, the general solution for the actual system 

Markov parameters y* = CA k B can be written as a 
convolution product of its previous values and the 
Markov parameters of the observer system as 


input single-output (MISO) problem. The results are then 
combined for MIMO identification. Multiple-input single- 
output problems, however, are essentially similar to 
single-input single-output (SISO) problems, therefore the 
algorithm can be best illustrated by first considering the 
SISO case. 


EXAMPLE 1: Consider a fourth-order single-input single- 
output system in observable canonical form driven by a 
random input sequence 



0 0 0 -0.305" 


r 

A = 

1 0 0 -0.110 

, B = 

0 


0 1 0 -0.110 


1 


-0 0 1 0.000- 


.1. 


C= Co 0 0 1] 


The first 30 Markov parameters are to be identified. The 
number of time step in each interval is thus chosen to be 
p = 30. With a 90-lime step input-output history, this 
results in 3 available repetitions for learning. In this 
example, the eigenvalues chosen are ±0.6 and ±0.7. In 
fact, any real, distinct value A, with magnitude less than 

one can be used as long as X? is sufficiently small to 
allow for accurate identification. In this example, the 
standard least squares method is employed, and 
convergence of the identified observer parameters is shown 
in Figure 1 . With these parameters, the observer Markov 
parameters are constructed, and then used to recover 
correctly the system Markov parameters. This is shown in 
Table 1. 


EXAMPLE 2: The algorithm is applied to identify a state 
space model of a mass-spring-dashpot system. This is a 
sixth-order system with three inputs and three outputs 
with the discrete-time system matrices given as 




0.9691 0.0154 
0.0154 0.9690 
0.0000 0.0077 
- 0.2817 0.1395 
0.1401 - 0.2838 
- 0.0005 0.0699 


0.0001 0.2120 
0.0155 0.0014 
0.9768 0.0000 
0.0009 0.9458 
0.1412 0.0180 
- 0.2122 0.0000 


0.0014 0.0000 
0.2139 0.0014 
0.0007 0.2127 
0.0180 0.0001 
0.9134 0.0182 
0.0091 0 . 9544 - 


Y k = CA k B 

= vl l) + X F/ 2 Vt. t ., 

i = 0 


( 45 ) 


Once the Markov parameters arc identified, a 
realization procedure such as ERA can be applied to 
obtain a realization of the system matrices. Physical 
aspects of the model such as natural frequencies, damping 
ratios, mode shapes can then be found. 


EXAMPLES 


In this section numerical examples are presented to 
illustrate the above developed identification algorithm. 
Recall that in this formulation, for the multivariable case, 
all the couplings between the outputs arc expressed in the 
"input” vectors. Hence for each output, the multiple-input 
multiple-output (MIMO) problem is treated as a multiple- 


B = 


0.0232 0.0001 0.0000 
0.0001 0.0233 0.0001 
0.0000 0.0000 0,0232 
0.2120 0.0014 0.0000 
0.0014 0.2139 0.0014 
- 0.0000 0.0007 0 . 2127 - 



0 0 0 0 0 

1 0 0 0 0 

0 1 0 0 0 . 


The first twenty five Markov parameters arc to be 
identified, thus p is chosen to be 25. Using a single time 
history of 100 time steps under random input excitation, 
this yields four repetitions for learning. Assuming for the 
moment, the true order of the system, n = 6 , is known. 
The eigenvalues for the observer are placed at ±0.2, 
±0.3, and ±0.4. The observer equation parameters 
corresponding to each output are first identified, and then 
the results combined to recover the actual system Markov 
parameters. Figure 2 shows the convergence of some 
arbitrarily selected observer parameter estimates 
corresponding to the first output. Results for the two 
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remaining outputs are similar, and not shown here. The 
algorithm correctly identifies all desired Markov 
parameters, each is a 3 x 3 matrix. 

The Eigensystem Realization Algorithm (ERA) is 
then applied to decompose the identified Markov 
parameters to obtain a set of realized system matrices 
which is equivalent to the set (A, B, C) above. They are 
listed here in that order as follows. 

0.9971 0.0561 0.0238 0.0606 - 0.0049 0 . 1166 ’ 

- 0.0553 0.9406 0.1879 - 0.1115 - 0.0291 0.0289 

- 0.0357 - 0.1762 0.9878 - 0.0341 0.1757 0.0753 

- 0.0497 0.1321 0.0173 0.9725 0.2365 - 0.0193 

0.0113 0.0232 - 0.1827 - 0.2027 0.9207 0.0923 

-- 0.1178 0.0050 - 0.0703 0.0326 - 0.1105 0 . 9600 . 

’ 0.2015 0.2770 0 . 2140 * 

- 0.1835 - 0.0138 0.3051 
0.0380 - 0.1269 0.2732 
0.2712 0.0430 - 0.1038 
0.1925 - 0.2085 0.0613 
L 0.1514 0.2515 0 . 1813 - 

0.1804 - 0.2296 0.2651 - 0.0293 - 0.1461 - 0 . 1934 ' 

0.2757 - 0.1481 - 0.1093 - 0.3034 0.1177 - 0.1246 

0.2309 0.2355 - 0.1408 0.1200 0.0625 - 0 . 2810 . 


Now, consider the identification algorithm used in 
Example 2, but this time the system order is assumed 
incorrectly. First, the syslcm is over-estimated to be of 
eighth-order. This calls for two additional eigenvalues in 
the algorithm, here chosen to be ±0.25. Note that since 
the assumed order of the system is higher, there are more 
parameters to be identified. However, at the final step, the 
algorithm correctly identifies the true order of the system 
and recovers all desired Markov parameters of the open 
loop syslcm. This means that over-parameterization docs 
not affect the final result. Second, the system is under- 
estimated to be only of second-order. This reduces the 
number of parameters to be identified substantially. Yet, 
again, the true order of the system and all the Markov 
parameters are recovered correctly at the final step. This 
indicates that under-estimation of the system order docs 
not lead automatically to incorrect identification. Correct 
results arc also obtained when the order of the system in 
this example is under-estimated to be 5, 4, or 3. The 
algorithm fails when the order is assumed to be 1, which 
is obviously an erroneous over-simplification. The 
identified observer parameters arc different under different 
order assumptions, yet they yield the same result in the 
final reconstruction of the actual system Markov 
parameters. 

EXAMPLE 3: In the absence of noise, the identified 
values arc found to be practically identical to the actual 
values. In the deterministic theory, this error can be made 
arbitrarily small simply by choosing p to be sufficiently 
large, and the eigenvalues A # sufficiently small such that 
the approximation in deriving Eq. (27) is valid. With 
regard to identification accuracy in the presence of noise, 
more effective recursive identification methods, e.g., the 


instrumental variable method, can be used to replace the 
standard least squares method at one step in this 
identification procedure. For example, consider die case in 
Example 1. When the output data is corrupted by about 
12% measurement noise, the standard recursive least 
squares method leads to biased results which can be 
corrected by using the instrumental variable method 
instead. Figure 3 shows this improvement in the accuracy 
of some arbitrarily selected identified observer parameter 
estimates. 

CONCLUSIONS 

This paper formulates an algorithm for identification 
of linear multivariable systems from a one set of input- 
output data. Instead of identifying the system directly, the 
proposed scheme identifies an observer for the system, 
whose poles can be assigned by an embedded real 
eigenvalue assignment procedure. Recursive techniques 
are used to estimate the matrix coefficients of an auto- 
regressive moving average model formed from this 
observer. These matrix coefficients are precisely its 
Markov parameters. From the identified observer Markov 
parameters, the true system Markov parameters are 
recovered through simple algebraic relations, and then 
used to obtain a realization of the system of interest. For 
modal identification, the modal parameters such as natural 
frequencies, modal damping, and mode shapes of the open 
loop system can then be found. Preliminary results 
indicates that the deterministic algorithm is fast and 
accurate. Stochastic aspects of the procedure is currently 
under investigation. 
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APPENDIX 


Table I: Identification of Markov Parameters of a Single-Input 
Single-Output System 


IDENTIFIED RECONSTRUCTED IDENTIFIED ACTUAL 


OBSERVER 

PARAMETERS 

OBSERVER 

MARKOV PARAMETERS 

MARKOV 

PARAMETERS 

MARKOV 

PARAMETERS 

-10.1021 

1.0000 

0.0000 



7,3337 

1.0000 

-0.9600 

-0.0005 


10.0710 

0.8499 

- 0.1100 

-0.0012 


-6.3025 

1.8500 

-0.9446 

-0.0015 

-0.0015 

-3.4628 


-0.0935 

0.0000 

0.0000 

2.6168 

1.3961 

-0.6336 

0.0023 

0.0023 

3.7143 

0.3142 

-0.0601 

0.0030 

0.0030 

-2.8682 



0.0031 

0.0031 


0.1707 


-0.0022 

imKi 


0.4850 

-0.2044 

-0.0080 

-0.0080 


0.0897 

-0.0188 

-0,0062 

-0.0062 


0.2605 

-0.1081 

-0.0051 

-0.0051 


0.0461 

-0.0099 

0.0114 

0.0114 


0.1359 

-0.0558 

0.0239 

0.0239 


0.0234 

-0.0051 

0.0077 

0.0077 


0.0695 

-0.0284 

0.0054 

0.0054 


0.0117 

-0.0026 

-0.0422 

-0.0422 


0.0351 

-0.0143 

-0.0652 

-0.0652 


0.0059 

-0.0013 

0.0135 

0.0135 


0.0176 

-0.0071 

0.0009 

0.0009 


0.0029 

-0.0006 

0.1331 

0.1331 


0.0088 

-0.0035 

0.1654 

0.1654 


0.0014 

-0,0003 

-0.1519 

-0.1520 



-0.0018 

-0.0079 

-0.0079 



-0.0002 

-0.3787 

-0.3787 



-0.0009 

-0.4029 

-0.4029 



-0.0001 

0.7800 

0.7800 


0.0011 

-0.0004 

-0.1100 

-0.1 101 


0.0002 

- 0.0000 

1.0000 

1.0000 


0.0005 

-0.0002 

1.0000 

1.0000 


1 o 








NUMBER OF TIME STEPS 


Figure 1: Convergence of identified observer parameters of a 

single-input single-output system 



NUMBER OF TIME STEPS 


Figure 2: Convergence of identified observer parameters of a 

multiple-input multiple-output system 



NUMBER OF TIME STEPS 


Figure 2: Convergence of observer parameters in the 

presence of noise using standard least-squares and 
instrumental variable methods 
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